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I. INTRODUCTION 

The recent world events have caused drastic changes in many areas and made the future 
that much more unpredictable. With this high degree of uncertainty and drastic financial 
cutbacks as evidenced in the Defense budget, the emphasis of technological advance and 
dominance is all the more heightened. With the requirement to "do more with less", it is 
imperative to find ways to save money. It is with this mind-set that better ways to prolong 
the life of operating machinery and defray unneeded maintenance costs be investigated. For 
many years, most transient machinery has gone unmonitored and removed at periodic time 
intervals for refurbishment regardless of the current machinery condition. Developing a 
method to monitor and diagnose these machines of a transient nature will provide a way to 
determine when overhaul or replacement is needed. 

Vibrational data has long been used as the tool to measure a machine's health. Early 
methods of monitoring machinery had been for the experienced personnel to listen or touch 
the machine to detect if a fault exists, however, this method is dependent on human senses 
and the same individual conducting the monitoring. It was for this reason that vibration levels 
as monitored by a dynamic signal-analyzer are now used for machinery health monitoring. 
The vibration levels of a good condition establish a baseline to allow a means of comparison 
when subsequent measurements are taken. These dynamic signal analyzers are time 
independent and do not show the time dependency required by transient machinery. It is 
because of this limitation, that a three-dimensional time, frequency and energy representation 
such as the Pseudo Wigner- Ville Distribution is being investigated as a method to monitor 


transient machinery. 


А. PREVIOUS WORK AND ACCOMPLISHMENTS 

This thesis is a follow-on work which has been preceded by LT G. Rossano and LCDR 
J. Calahorrano, both in conjunction with Prof. Y.S. Shin. Rossano established the initial form 
of the computer program that calculates the discretized Pseudo Wigner-Ville Distribution. 
He evaluated the effect that different signals and procedures had on the Distribution along 
with characterizing an actual pump signal. Calahorrano made some further improvements to 
the computer program and then applied this improved version of the computer program to 
characterize some classified signals. Both of these works have been a great benefit in dealing 
with characterizing signals, but have not investigated the smoothing method and energy or 


amplitude of the discretized Pseudo Wigner- Ville Distribution. 


B. OBJECTIVES OF CURRENT RESEARCH 

The Pseudo Wigner- Ville Distribution provides an outstanding time-frequency domain 
spectra representation of an input signal as shown through the works among others of Rossano 
and Calahorrano here at the Naval Postgraduate School. The purpose of this research is: 


• Toreview and modify the computer program developed here at the Naval Postgraduate 
School to support the smoothing window and energy analysis work. 


e To conduct a smoothing window sensitivity analysis. 
• To investigate the energy content of the discretized Pseudo Wigner- Ville Distribution. 


e To apply the developed smoothing window and energy principles to an actual 
machinery monitoring application. 


e To investigate the feasibility of using the Pseudo Wigner- Ville Distribution energy as 
a pre-processed input to a Neural Network. 


Il. THE PSEUDO WIGNER- VILLE DISTRIBUTION 


A. THEORETICAL BACKGROUND 
1. Evolution 

The current day Pseudo Wigner-Ville Distribution is a three dimensional (time, 
frequency, and energy) representation of a signal that is particularly well suited for analysis 
of non-stationary signals. Eugene Wigner [Ref. 1] first introduced the Wigner distribution in 
1932 to study the problem of statistical equilibrium in the area of quantum mechanics. His 
work was furthered in 1948 by J. Ville [Ref. 2], who used the Wigner distribution in the area 
of signal analysis. A major contribution of Ville’s, is the use of analytic signals as an input 
to the distribution vice the customary real signal. The advantages of using an analytic signal 
in the distribution are two-fold [Ref. 3]. First, the distribution of a real signal results in a 
symmetrical spectrum with only half of the result containing useful information, while the 
use of an analytic signal avoids this negative frequency redundancy. Secondly, by accounting 
for only the positive frequencies it satisfies both the practical and the mathematical 
completeness of the problem. The third part of the Pseudo Wigner-Ville Distribution (ie. 
Pseudo) , arises from the need to apply a smoothing window to the resulting distribution. 
This is a result of the presence of cross terms that arise from multiple component signals and 
will be covered in depth later. 

The works of T.A.C.M. Claasen and W.F.G. Mecklenbrauker [Refs. 4,5, and 6] in 
1980 have paved the way for many recent applications of the Pseudo Wigner- Ville 
Distribution. Their three part paper is an all encompassing work that has served to highlight 
the capabilities of the Pseudo Wigner-Ville Distribution and allow for greater knowledge of 


the subject with an emphasis on signal processing and analysis. 


A considerable amount of recent work using the Pseudo Wigner- Ville Distribution 
have come in the areas of optics [Refs. 7,8, and 9] and speech [Ref. 10 and 11]. Additionally, 
T.J. Wahl and J.S. Bolton of Purdue University [Ref. 12] have used the Pseudo Wigner- Ville 
Distribution to analyze structural impulse responses. Two recent papers have investigated the 
use of the Wigner- Ville Distribution in the area of machinery monitoring. Flandrin et. al. 
[Ref. 13] proposed the Wigner-Ville as a means to confirm a machinery diagnosis in the 
specific application of a Pressurized Water Reactor power plant and lastly, Forrester [Ref. 14] 
has investigated the Wigner-Ville Distribution as a method for fault detection in gears. As 
evidenced by the cited examples, the capability and adaptability of the Wigner- Ville 
Distribution is enormous with applications of it's use rapidly expanding. 

2. Function Definition 
The Wigner-Ville Distribution is a transformation of a signal into the time- 


frequency domain. By definition, the cross Wigner- Ville Distribution is defined as: 


WDF, s= |e r(t+=) s*(t->) dr 


where: r(t) = a complex time history 
s(t) = a complex time history 
t = time 
w = frequency 
* = complex conjugate 


Similarly, the auto Wigner-Ville Distribution is defined as: 


WDE, La, Jeer riëa zl к'(е-—) ф 


Since the concentration of this research is with the representation of a single input signal, the 
auto Wigner- Ville Distribution will be used. This expression is essentially a Fourier transform 
of the auto-correlation of a signal, which may be thought of as a three dimensional spectral 
density function. 
3. Distribution Properties 

There are several properties of the Distribution that are worth noting. The 
properties listed below have been shown in many works concerning the Wigner-Ville 
Distribution. These properties will be shown in the form of the auto Wigner-Ville 
Distribution [Ref. 15]. 

Time and frequency shift properties follow: 


Ф А time shift in the signal corresponds to a time shift of the Distribution: 


DEE (t, ©) = ИРЕ "ашиг ©) 


Ф A frequency shift in the signal corresponds to a frequency shift in the Distribution: 


WDF 3025 засе (Ё, о) = WDE: (Е, (у - 03 ) 


e Both a frequency and time shift in the signal corresponds to the same in the 
Distribution: 


ИН ос 22) 28:22) Ce о) ын WDE: (CET, о-0) 


Since our concern is in the area of machinery monitoring and diagnostics, the aforementioned 
Ner 2 ZE г 4-7 


properties are a nessesity in the development of a suitable monitoring program. 


The next three properties provide information concerning the energy contained 
within the Wigner- Ville Distribution. 


• The integration of the Distribution over the frequency domain provides the 
instantaneous signal power: 


1 
за ЈУРЕ. (е, а) ае -| s(t) P 


e The integration of the Distribution over the time domain provides the energy density 
spectrum: 


| кре, (с,ө) 46 -18(о) Р 


- 00 


e The integration over both the time and frequency domain provides the total energy of 
the input signal: 


в 
= | | пре, „„(5, в) асао - 11511 


- 00- со 


Again, as with the time and frequency properties, the energy contained in the Wigner- Ville 
Distribution must be definable in order to use this in a machinery monitoring program. 
Additional properties may be found in Reference 4. 
4. Cross Terms Within The Distribution 
A particularly troublesome characteristic of the Wigner-Ville Distribution is the 
presence of cross terms. These cross terms result in the Distribution taking on negative values 
as shown later in Figure 2. Since the properties above have shown that the Wigner- Ville 


Distribution represents an energy value, these negative values raise an area of concern. Jones 


and Parks [Ref. 16] show that these cross terms are essentially the interference of multiple 
component signals that occur when there is a non-zero cross-correlation within the signal. 
The equation listed below shows the cross-term that results from a signal that consists of two 


components. 


WDF,,,(€,@) - WDF,(t,@) + 2RE[WDF, ,(t,@)] + ЮРЕ, (Е, о) 


х+у 


As shown, the cross-term is represented by the center term while the other tems represent 
their respective signal components. There are several methods for eliminating these terms and 
providing a more presentable Distribution and there has been a considerable amount of 
research in the are of an optimum smoothing method for the Wigner-Ville Distribution as 
indicated by References 3,4,17, and 18. The method chosen for smoothing the Wigner- Ville 
Distribution in this work is a post processing, sliding Gaussian window function. What is 
meant by post processing, is that the distribution is smoothed after calculation of the Wigner- 
Ville Distribution. 
5. Gaussian Smoothing Function 

The Gaussian window function was first proposed mathematically as a means to 

smooth the Wigner- Ville Distribution by N.D. Cartwright [Ref. 19]. The continuous form of 


the Gaussian window function follows: 








G(t,@) = ———e 


Wahl and Bolton in Reference 12 have furtherd this work and developed the sampled form 


of the window function that is used in our discretized form of the Pseudo Wigner- Ville 


Distribution for smoothing. The effect of various size smoothing windows will be presented 


in greater depth in Chapter III. 


В. WIGNER-VILLE DISTRIBUTION COMPUTER PROGRAM 
The computer program listed in the Appendix is a method that was developed here at 
the Naval Postgraduate School for calculating the discretized Pseudo Wigner- Ville 
Distribution. As mentioned earlier, the program that computes the discretized Pseudo 
Wigner-Ville Distribution is largely the work of G. Rossano [Ref. 15], who initially 
investigated the use of the Wigner- Ville Distribution as a tool for machinery monitoring. For 
the purposes of this work, the program has been thoroughly reviewed and rewritten to allow 
for ease of readibility and modification. The current version enables the user to look at not 
only the resulting smoothed version of the Pseudo Distribution, but also the unsmoothed 
Wigner-Ville Distribution. The second major modification was to program the smoothing 
window subroutine to allow for the use of different size windows in smoothing the 
Distribution. These changes have proven invaluable for evaluating the effect of various sized 
smoothing windows and in conducting the energy sensitivity analysis. Presented below are 
the two major equations that comprise the calculation of the discretized Pseudo Wigner- Ville 
Distribution as performed in the program listed in the Appendix. 
1. Discrete Time Wigner-Ville Distribution 
As mentioned earlier, the Wigner-Ville Distribution is essentially a Fourier 
transform of an auto-correlation function. There have been a couple of equations developed 
to represent the discrete Wigner-Ville Distribution. The version that has been used in the 
development of the program listed in the Appendix, is one that was developed in Reference 
4 by Classen and Mecklenbrauker and and later used in the work of Wahl and Bolton [Ref. 
12] and Rossano [Ref. 15]. The only difference of note between the various discrete versions 


is essentially a difference in constants. All of the versions are basically a Fourier transform 


of an auto-correlation function. Shown below is the discretized equation that has been used 


for programming purposes in this work: 


N 


WDF, ,(1à9,jAt) - Y, RelFFT[2A t СОКК (1) ]) 
7-1 


where: Re = Real part 
Corr(i) = complex auto-correlation of the signal s(t) 
FFT = Fast Fourier Transform 


At = sample time 


2.  Discretized Gaussian Window Function 
Due to the need to deemphasis the resulting cross-terms within the Wigner- Ville 
Distribution, a suitable method for smoothing has been obtained. The sampled equation of 


the selected smoothing window is shown below [Ref. 12]: 


21 (тас)? , (лав)? 
Ж 2 2 
Gt 28 о») e Боа: =: ме 20. 292 
270 ,0, 
where: -2М<т<2М & -2N<n<2N 


о, = ММ & os NA 


Ш. SMOOTHING WINDOW ANALYSIS 


As shown above, the smoothing function or window that is used in this work is a 
Gaussian window that is applied to the unsmoothed Wigner- Ville Distribution. Figure 1 15 (ће 
raw time input signal of a 100 & 400 Hz sine wave that will be processed through the Wigner- 
Ville Distribution. Figure 2 1s the unsmoothed Distribution that results from this combined 
100 and 400 Hz sine wave. As can easily be seen in Figure 2, the resulting cross terms are 
very dominant and need to be deemphasized. This Gaussian smoothing function is applied 
after the Wigner- Ville Distribution has been calculated (1e. post processing), which allows for 
visual inspection of both the unsmoothed Wigner-Ville Distribution and the effect that 
different sized windows may produce on the smoothed Distribution (ie. Pseudo). For 
example, Figure 3 is the smoothed result of Figure 2 after applying a 13 by 13 sized window. 
The next few sections will discuss some of the characteristics of this post-processing 
smoothing window and what may be gained or lost through varying the size of the applied 


window. 
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Figure 1. 100 & 400 Hz Sine Wave Input 
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Figure 2. Unsmoothed Wigner- Ville Distribution 
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А. SIZE CRITERION 

As shown by N.D. Cartwright [Ref. 19] for the continuous form of the Wigner- Ville 
Distribution, the size of the smoothing window may be selected such that a positive 
Distribution will result. Following, are these criterion that have been mathematically proven 


for the continuous case to assure a positive distribution: 


2 


1 
а 0 > — 
і 4 


en 


where: 


с, = МАГ с, = NAO 


By chosing M and N such that this criterion is met, the resulting smoothed discretized 
Wigner- Ville Distribution does not necessarily result in all of the values being positive. This 
may be seen in Figure 4, where the window that was used was again a 13 Бу 13 (іе. М = 13 
& M = 13) sized window. However, it does deemphasize the large cross-terms to the point 
that they are nearly zero and yields a very presentable Pseudo Wigner- Ville Distribution. The 
ability to deemphasize these cross terms is a necessity if the desire is to evalute the energy 
change within a signal. For the remainder of this work, the window size selected (ie. N & M, 


N X M, or N by M) will be such that this criterion is met as closely as possible. 
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Figure 4. Smoothed Distribution with Negative Values 


В. TIME AND FREQUENCY DEPENDANCE 

Since the criterion shown previously are time and frequency dependent, there must be 
a relation between the two that can be reduced to one parameter. The customary time and 
frequency relations of signal analysis are slightly changed in the context of the Wigner- Ville 
Distribution. Shown below are the time and frequency relations when working with the 


Wigner- Ville Distribution. 


, * DPxAC ; Gm Жж ; - 2xDPxAf 


t Е 


where: пах = Maximum signallength 4 бах = maximum frequency 


m m 


At = sampling time & Af = frequency resolution 


DP = number of data points 


The definitions of the size criterion discussed earlier for the smoothing window follow: 


202 L 
С t9o 2 4 
and since: 
02 - (МАБ)? o = (NAw)? 


and the time and frequency are related in the Wigner- Ville Distribution through the number 


of data points as shown: 


=^ 2 DP At 


Making this substitution into the size criterion, yields a very simple result for determining the 


size of the smoothing window required. 


MxN 2 DP 
л 


In looking at this result it seems as if the window is completely independent of time and 
frequency, but it must be rememberd that M and N are parameters that specify the lengths 


of the smoothing window in the time and frequency directions respectively. 


C. RELATION TO NUMBER OF DATA POINTS 

Another area of interest is the relation of the window size to the number of data points 
in the sampled signal. The best way to show this relation is through an example. For 
instance, a signal of a length 4.0 seconds that is sampled using 512 data points, will result in 
a At = 7.8125 msec and Af = 0.0625 Hz as calculated using the constraints of time and 
frequency as defined by the Wigner- Ville Distribution. Since the example has used 512 data 
points, the previous section has shown that chosing N = 10 and M = 16 will approximately 
satisfy the size criterion limitations. Additionaly, by definition of the smoothing window, the 


area of coverage will be: 


-2М<т<2М -2N<n<2N 


This means that for the example shown, the coverage of the window will be 40 steps in the 


frequency domain and 64 steps in the time domain. 


40хАГ - 2.5 64xAt = 0.5 AREA = 1.25 


where: АЁ =0.0625Hz & А = 7.8125 тѕес 


The AREA = 1.25 is a unitless quantity that defines the area that the smoothing window 
encompasses while smoothing one point in the Distribution. The unitless dimension is the 
result of the time being in seconds and the frequency in hertz. For comparison purposes, the 
total area of the distribution is the product of f nax aNd tmax Which is equal to 256.0, resulting 
in the window smoothing over 0.488% of the total Distribution, while smoothing one point. 

Now consider another example that uses 1024 data points. As shown earlier, the 
product of M and N must be greater than or equal to the number of data points divided by 
pi in order to meet the window size criterion as established by Reference 12. With respect to 
the earlier example, the size of the window (ie. N times M) must increase two times. 
Therefore, a selection of N = 10 and M = 32 will be used for this example. These values are 
very Close to the required criterion and will suffice for this illustrative example. The signal 
length (tnax) Will be 2.0 seconds resulting in fi, = 256.25 Hz, At = 1.953 msec, and Af = 


0.125 Hz. In relation to the area of coverage of the window, the following results apply: 


40xAf = 5.0 128xAt - 0.249984 AREA « 1.25 


where: Af =0.125 Hz & At = 1.953 msec 


This is the same result as the 512 data point case, with the following exception. The total area 


of the Distribution (tnax X f, a x = 512.0) has doubled, while the area of the smoothing window 


m 
coverage has remained the same. This results in the window only covering 0.24496 of the 
entire Distribution when smoothing a single point within the Distribution. In conclusion, the 


greater the number of sampled data points within the processed raw time signal, the smaller 


the area the smoothing window will encompass while smoothing the Distribution. An 
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important trade-off to consider, however, is the increase in computer time required to run 


the additional number of data points. 


р. TIME/FREQUENCY RESOLUTION 

Since the smoothing window used in this work gives the ability to vary the time and 
frequency size of the window, the following examples show the results that can be obtained. 
Figure 5, is the input signal that is sampled using 512 data points. It is a 400 Hz sine wave 
that has been computer generated with an amplitude of 2.0. In order to be consistent with the 
previous work, the size criterion established earlier will be used as a guide for the selection 
of M and N in the smoothing window. Figure 6, is the result of a window in which M and 
N are both equal to 13. Since this smoothing window does not require M and N to be equal, 
varying these parameters yielded some intersting results. Figure 7, shows the result of the 
input signal being smoothed with a 35 by 5 window. This means that N = 35 and M = 5, or 
in other words, the window is longer in the frequency domain than it is in the time domain. 
The resulting distribution has an amplitude that is much more constant throughout the entire 
time length, but the amplitude itself has been reduced considerably as compared to the 13 by 
13 example. Figure 8 is just a different view of this example, showing the loss of frequency 
resolution with a 35 by 5 window. Just the opposite case to this example is where the window 
used for smoothing is a 5 by 35 window. In this case, the window is longer in the time 
domain than the frequency domain and just the opposite results are obtained. Figure 9, shows 
the result of this window, while Figure 10 gives a view showing the degree of frequency 
resolution that can be obtained. This can be very beneficial in the monitoring of machinery, 
since in many cases the ability to differentiate between signals that are very close to one 


another is paramount. 
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Figure 5. 400 Hz Sine Wave Input 
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Figure 7. Pseudo Wigner-Ville Distribution (35 X 5 Window) 
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Figure 8. Frequency Resolution (35 X 5 Window) 
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Figure 9. Pseudo Wigner- Ville Distribution (5 X 35 Window) 
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Figure 10. Frequency Resolution (5 X 35 Window) 
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E. TRANSITIONAL PEAK PHENOMENON 

In works that have given examples of the Wigner-Ville Distribution for transient 
signals, the point at which the signal’s frequency is changing from increasing to decreasing 
or vice versa has a large amplitude increase in the Wigner- Ville Distribution. For example, 
Rossano [Ref. 15] shows an artificially generated linear chirp and also an actual pump 
signature that has what appears to be a large concentration of energy at the transition point. 
In order to use the Wigner- Ville Distribution for transient machinery monitoring, there must 
be an explanation or way to quantify this "ghost peak" that 1s occuring at the transition point. 
To date, there has not been an explanation to this phenomenon, however, by looking at the 
unsmoothed Wigner- Ville Distribution, the reason for this is easily seen. Figure 11 is a linear 
chirp sine wave that is used as the input to the Wigner- Ville Distribution. Figure 12, shows 
the resulting unsmoothed Wigner- Ville Distribution of this linear chirp and as discussed and 
shown earlier, the cross-terms can be seen occuring exactly at the mid-point of the conflicting 
frequencies and times. These cross-terms build until they reach a maximum concentration 
at the transition or frequency change point. Since these are nothing more than cross-terms, 
there must be a way to deemphasis these through the judicial use of the smoothing window. 
Figures 13 through 16 show the result of varying the smoothing window from a point of 


essentially no effect to entire elimination of the ghost peak. 
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Figure 11. Linear Chirp Sine Wave Input 
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Figure 12. Unsmoothed Wigner- Ville Distribution 
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Figure 13. Pseudo Wigner- Ville Distribution (5 X 35 Window) 
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Figure 14. Pseudo Wigner- Ville Distribution (13 X 13 Window) 
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Figure 15. Pseudo Wigner- Ville Distribution (18 X 10 Window) 
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Figure 16. Pseudo Wigner-Ville Distribution (35 X 5 Window) 
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This chapter has shown the benefits that may be gained through the use of the 
smoothing window. The frequency resolution and amplitude are the two major areas affected 
by the window changes. The last area presented, showed the reason for the large amplitude 
increase at the frequency transition point and how these "ghost peaks" may be smoothed and 
to what extent deemphasized by adjusting the smoothing window size. Within the next 
chapter, a thorough energy analysis of the discretized Pseudo Wigner- Ville Distribution will 


be conducted. 
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IV. WIGNER-VILLE DISTRIBUTION ENERGY ANALYSIS 


А. | ENERGY DEFINITION 

As mentioned earlier, the Wigner- Ville Distribution is essentially the Fourier transform 
of the auto-correlation of an input signal. This may be thought of as the spectral density 
function in a three dimensional sense. The total energy of the Wigner- Ville Distribution was 


listed earlier as: 


-1_[[ирр(о, t) dwdt - y? 


= 09-- œ 


where: 


y? - mean square value of the signal 


There are a few points that must be made concerning this equation. The first is that the 
integration is from minus infinity to plus infinty and another is that it is the integration of 
the resulting Distribution of a real signal, not of an analytic signal as used 1n most discretized 
forms of the Wigner-Ville Distribution. Lastly is that this equation has been proven 
mathematically, which unfortunately does not take into account the fact that there are 
resulting cross-terms and that the Distribution is subsequently smoothed because of this. 
Since the last point is indicative of what is seen in actual real world problems, the energy must 
be defined in a discretized sense that accounts for the above mentioned points and also these 


unavoidable cross-terms. 
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B. EFFECT OF WINDOW SIZE ON ENERGY 

Since smoothing of the Distribution is a necessary requirement, the size of the applied 
window is something that must be decided upon. Another important point is knowing what 
effect this will have on the energy of the Distribution. In the previous chapter, the size 
criterion that was established by Cartwright [Ref. 19] and later used by Wahl and Bolton in 
Reference 12 was discussed. Within this discussion, the visual benefits of varying the window 
size in the time and frequency domain were presented. In order to define the energy, the 
effect on the energy of the Pseudo Wigner-Ville Distribution as caused by varying these 
windows will be presented. 

Figure 17 1s a 400 Hz sine wave of length 0.256 seconds and amplitude 2.0 that will be 
sampled and have the Wigner-Ville Distribution calculated. Figure 18 is the resulting 
smoothed Wigner- Ville Distribution that has been processed using a 13 by 13 sized smoothing 
window. Figure 18 also shows that the volume is equal to 0.0709118 V?. This is 
representative of the energy contained within the Wigner- Ville Distribution of the signal in 
Figure 17 and is nothing more than the volume under the resulting Distribution. Ву 
definition, this volume should be equal to the mean square value of the 400 Hz sine wave, 
however, due to the points mentioned above, this is not going to hold true and will be 
discussed in depth later. Figures 19 and 20 are also smoothed Distributions of this same 400 
Hz signal only using a 35 by 5 window in the case of Figure 19, and a 5 by 35 window in 
Figure 20. Since N times M (the size criterion) are very close in these three examples, the 
resulting energy values of the three cases are essentially the same. The point to be emphasized 
from this is that as long as the area of the window does not vary significantly, the energy will 
be represented in a constant manner and independent of the time and frequency window size 
selected. The last example that will be presented in this section is that this result 1s 
independent of the frequency of the signal as shown by the combined 100 and 400 Hz sine 


wave example of Figure 21, in which the energy is equal to 0.1418314 V?. This is essentially 
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twice the value of the single sine waves of Figures 18 through 20. Of special importance for 
transient machinery applications is Figure 22, which is a transient or ramped sine wave that 
shows the same resulting energy as the previous stationary signals. This shows that the 
Wigner-Ville Distribution energy is independent of the nature of the signal allowing for the 
possibility of monitoring a piece of machinery during both it’s steady state and transient 


phases of operation. 
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Figure 17. 400 Hz Sine Wave Input 
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18. Pseudo Wigner-Ville Distribution Energy (13 X 13 Window) 
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Figure 19. Pseudo Wigner- Ville Distribution Energy (35 X 5 Window) 
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Figure 20. Pseudo Wigner- Ville Distribution Energy (5 X 35 Window) 
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Figure 21. Combined Sine Wave Distribution Energy (13 X 13 Window) 
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Figure 22. Ramped Sine Wave Distribution Energy (13 X 13 Window) 
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The previous figures show that if the Wigner- Ville Distribution is going to be used to 
accurately portray or maintain signal information, then the window size must at least remain 
constant and better yet to keep the identical window size for each measurement of 
comparison. This requirement to maintain a constant sized smoothing window is not the only 
constraint that must be meet for energy analysis purposes. The next section will discuss 


another of these elements. 


C. TIME DEPENDENCE OF WIGNER-VILLE ENERGY 

Not only is the window size going to influence the energy result, but the time length 
of the signal will also play a large part. This varies considerably from the conventional 
knowledge of the mean square value of a signal being independent of time. The standard 
mean square value will be the same for a signal of time length 2.0 seconds as it will for a 
signal of 4.0 seconds. For instance, the mean square value of a sine wave x(t) is: 


1 2 
y? - [ (auto spectrum) df -= 2- 
0 


where: x(t) = X sin (ot) 


To compute the mean square value of this sine wave, it is independent of time and will equal 
0.5 Гог Х = 1.0. This is not the case for the mean square value or total energy of the signal 
as found using the discretized form of the Wigner-Ville Distribution. The following figures 
will show this second important point that must be realized when comparing Wigner- V ille 
energies of various signals. Figure 23 is the 400 Hz input signal of length 0.128 seconds. This 
time length is half of that of the 400 Hz signal used in the previous section. The resulting 
Wigner-Ville Distribution is shown in Figure 24 with an energy value of 0.0354558 А.С 


can easily be seen, this resulting energy value is not the same and in fact, the energy of Figure 
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24, is one half that of the earlier 400 Hz signals of length 0.256 seconds. Had the signal been 
of a length 0.512 seconds or twice as long as the earlier examples, the resulting Wigner- Ville 
energy would have doubled. The window size used for this Pseudo Wigner- Ville Distribution 


was a 13 by 13 sized window. 
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Figure 23. 400 Hz Sine Wave (Tj, - 0.128 sec) 
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This has yielded a second concern that must be kept in mind when using the energy of 
the discretized Pseudo Wigner-Ville Distribution. Since the energy is proportional to the 
signal length, an easy conversion can be made by using the lengths of the individual raw input 
signals. This time dependence holds true for not only steady state signals as shown, but also 
for signals of a transient nature such as previously shown with Figure 22. Throughout these 
last two sections, the comparison of energy values of signals has been mentioned a great deal. 
The thought process behind this is that if there is an amplitude change in the input signal, this 
will be reflected in the resulting Pseudo Wigner- Ville Distribution. By ensuring that these 
previous two conditions are met, an investigation into the resulting energy changes from an 


amplitude increase or decrease of an input signal can be made. 


D. ENERGY PROPORTIONALITY OF THE DISTRIBUTION 

The mean square value of a simple sine wave is the amplitude squared divided by two. 
This was shown in an earlier equation and now will be used to show the proportionality of the 
Wigner-Ville Distribution. For instance, the following applies for signals of differing 


amplitudes as shown: 


amplitude - 1.0 ў? “ 0.5 


amplitude “ 2.0 Ш^ = 2.0 


This amplitude increase from 1.0 to 2.0 represents an increase of four in the mean square 
value of the signal. Figures 25 and 26 are the raw input sine waves of amplitude 1.0 and 2.0 
respectively. The resulting Pseudo Wigner- Ville Distributions are shown in Figures 27 and 
28, which acurately portray the four-fold energy increase. These Distributions have been 


smoothed using a window of 13 by 13. 
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Figure 25. 400 Hz Sine Wave (Amplitude = 1.0) 
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Figure 26. 400 Hz Sine Wave (Amplitude = 2.0) 
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Pseudo Wigner-Ville Distribution (Amplitude = 1.0) 
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ure 28. Pseudo Wigner-Ville Distribution (Amplitude = 2.0) 
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These previous four figures show the proportionality of the energy within the Wigner- 
Ville Distribution. These two cases are just an example to show the proportionality of the 
energy. Runs were made for many more cases (ie. X = 4.0, X = 8.0, etc) and the 
proportionality relation held true as shown with the previous Figures. This was another 
required element in order to use the energy as the basis for a method of machinery 
monitoring. Had the cross-terms increased in such a fashion as to not allow the Distribution 
to show the proportional energy change, the energy content of the Pseudo Wigner- Ville 
Distribution would have been rendered useless for the purposes of comparison and evaluation. 
The next chapter will use these requirements and basics that have been developed to show an 
actual application of how the Wigner-Ville energy may be used to conduct machinery 


monitoring. 
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У. GEAR MODEL APPLICATION 


A. DESCRIPTION OF GEAR MODEL 

The simple complexity gear model used in this work was built based upon a model used 
by D.K. Carlson [Ref. 20] for Artificial Neural Network applications for machinery 
monitoring. It may be of a simple nature, however, it provides all of the necessary 
components that may be found in more complex systems. A schematic of the gear model can 
be seen in Figure 29. The two major changes of this model over the prototype used by 
Carlson is that the components are much larger enabling a better signal to noise ratio and 
secondly, the motor used to drive the gear model 15 mounted on a separate base to reduce 
motor enduced vibrations. 

The model consists of a single reduction gear train, with the pinion gear being a 90 
tooth spur gear and the driven gear a 120 tooth spur gear. The gears were manufactured by 
Martin Corporation with a 1/2 inch bore and a 14.5 degree pressure angle. Aluminum blocks 
housed the NICE 1/2 inch bore radial ball bearings that were used to support the shaft and 
gears. These aluminum blocks were bolted to the support base which was a 1.0 inch thick 
plexiglass slab. Driving the single shaft is a General Electric 1/6'^ HP motor, that is 
controlled by the use of a Bodine Electric Company combination rectifier and variable 
potentiometer speed controller. The shaft speed was monitored through use of a Power 


Instruments Model 1720 RPM Optical Proximeter. [Ref 20] 


B. SIGNAL SAMPLING EQUIPMENT 
This section will describe the components shown in Figure 30, that were used to sample 
the vibrational signal of the gear model. The signal was monitored using an ENDEVCO 


model 303A03 accelerometer that was amplified and powered with a PCB model 483B07 
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power unit. The signal was then processed through a Krohn-Hite model 3342 analog filter 
configured in a band pass mode, prior to being sent to the HP 3565S system for sampling and 
Storage. Lastly, the signal was transfered to the VAX 3520 workstation for the calculation and 
plotting of the Pseudo Wigner- Ville Distribution. 

1, ENDEVCO Model 303A03 Accelerometer 

The ENDEVCO Isotron PE accelerometer is a minature accelerometer that has a 
medium range high frequency capability. The operation of the accelerometer is based upon 
a piezoelectric quartz transducer sensing element. The sensitivity of the accelerometer is 10 
mV/g with a resonant frequency of 70 kHz. The resolution available is 0.02 g, with a 
maximum range of + 500 g. 

The accelerometer was mounted on top of the aluminum bearing housing that 
supports the shaft and specifically at the position shown in Figure 30. The accelerometer wass 
attached by means of a mounting wax which yields a maximum operating range out to 15k - 
20k Hz. This accelerometer signal was then sent through the power unit described below. 

2. PCB Model 483307 Power Unit 

This power unit is used to power the low impedance quartz transducers and 
amplify the signal if desired. The power unit provides an adjustable 2 to 20 mA current for 
purposes of transducer excitation. The gain adjustment is available from 0 to 100 and is set 
through the use of a ten turn vernier gain pot and a three position gain multiplier switch. For 
the purposes of these tests, the gain was set to 20 before sending the signal to the analog filter. 

3. Krohn-Hite Model 3342 Anolog Filter 

The model 3342 variable filter is a digitally tunned filter that will function as a 
High-Pass or Low-Pass filter. When the two channels of the filter are connected together, the 
filter will function as a band pass filter, which is the configuration for the gear model work. 


The range of the filter is from 0.001 Hz to 99.9 kHz as adjusted by three rotary decade 
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switches and a rotary six positiohn multiplier switch. In addition, the filter unit has a gain 
setting of unity (0 dB) and 10 (20 dB). The signal was not further amplified using this 
capability of the filter prior to sampkling by the HP 3565S system. 
4. НР 3565S Measurement Hardware System 

The HP 3565S measurement hardware with the HP- Vista software uses a HP-9000 
series computer system for controlling purposes. The HP 3565S system is a modular 
multichannel system that can process data in both the time and frequency domain. The 
system is capable of handling 64 source and input modules, but is presently configured for 
eight. The HP-Vista software allowed for the sampling and storage of the preprocessed 
accelerometer signal that was later processed through the Pseudo Wigner- Ville 


distribution. 
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Figure 30. Signal Sampling Equipment Schematic 
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с. ESTABLISHING MONITORING FREQUENCIES 

This section will describe how the monitoring frequencies have been determined and 
how they will be used to monitor the gear model and in future applications a piece of 
transient machinery. As shown earlier in Chapter III, the energy contained in a signal is 
accurately portrayed within the resulting time-frequency domain of the Wigner-Ville 
Distribution. Through the Wigner-Ville Distribution energy that is contained in the 
established monitoring frequencies, faults may be detected. The shaft frequency is just the 
speed of the shaft, while the gear mesh frequency is the number of teeth on the rotating gear 
multiplied by the shaft speed or frequency. Both of these frequencies have associated 
harmonics that can be of equal importance as the frequencies themselves in the monitoring 
of machinery, It is for these reasons that the frequencies listed below have been chosen based 
upon a shaft frequency of 15 Hz. 

1. Monitoring Frequencies 

For the purposes of this work, just the shaft and gear mesh frequencies of the 

model will be included. Just as in standard established machinery monitoring methods, the 
shaft frequency and associated harmonics can provide a great deal of information for fault 
detection and diagnosis. Along with the shaft frequencies, the gear meshing frequencies 
provide valuable information of faults and will also be monitored. The purpose of the two 
examples will be to show that the Wigner-Ville Distribution energy (ie. volume) changes as 
a result of damage and will prove to be a very valuable tool in future machinery monitoring 
methods. The established monitoring frequencies have been filtered as follows: 

e Shaft Input 1: (5 - 100 Hz) This is the spectrum of signals in the lower frequency range 
selected to determine the total energy contained within the shaft frequency and it's 
harmonics. 

e Shaft Input 2: (10 - 20 Hz) This frequency range captures just the shaft frequency. 


+ Shaft Input 3: (25 - 35 Hz) This frequency range captures the 1** shaft harmonic. 
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+ Shaft Input 4: (40 - 50 Hz) This frequency range captures (ће 279 shaft harmonic. 


e Shaft Input 5: (55 - 100 Hz) This frequency range captures the upper harmonics of the 
shaft frequency. 


e Gear Mesh Input I: (1250 - 1450 Hz) This frequency range captures the gear meshing 
frequency. 


e Gear Mesh Input 2: (2600 - 2800 Hz) This frequency range captures the 1** gear mesh 
harmonic. 


« Gear Mesh Input 3: (3950 - 4150 Hz) This frequency range captures the 2"¢ gear mesh 
harmonic. 


2. Combined Monitoring Technique 

Using the Wigner- Ville energy as an input to an Artificial Neural Network is a 
follow-on goal of this research. It is for this reason, that the above established frequencies 
have been identified as inputs. Within this scheme, the Wigner- Ville is used to characterize 
the signal and give a single energy value at the particular frequency of interest as an input to 
an Artificial Neural Network. By filtering these signals prior to processing, just the 
frequency of interest will be captured along with it’s respective Wigner- Ville energy content. 
The Wigner- Ville energy data can now be used to produce a comprehensive training set for 
the Artificial Neural Network. By pre-processing the inputs to a Neural Network with the 
filtering and the Wigner-Ville Distribution, the number of inputs is greatly reduced. Since 
these associated energy levels are fairly stable for a given machinery condition, the Neural 
Network will be capable of detecting patterns in energy level shifts that are associated with 
particular machinery faults. To support this last statement, the following section will show 
two examples of machinery faults imposed on the gear model and the response of the Wigner- 


Ville energy to these faults. 
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D. MACHINERY FAULT EXAMPLES 
1. Mechanical Looseness 

This fault was imposed on the gear model by loosening the aluminum support 
blocks that house the bearings and support the shaft. Specifically, the blocks that support the 
driver shaft were loosened. The degree of looseness imposed was difficult to judge and made 
this particular fault a difficult one to diagnose. The blocks were not completely undone, but 
loosened just a small amount from their original firmly bolted posture. Shown in the table 
below are the results of averaging two good runs and two post-damage runs and comparing 
the respective energy values at the different frequencies. These energy values are nothing 
more than the volume under the Pseudo Wigner- Ville Distribution as discussed earlier. Only 
the shaft frequencies are shown for purposes of clarity and since it is in these monitored 


frequencies that a fault of mechanical looseness should present itself. 
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Table I. MECHANICAL LOOSENESS 


MECHANICAL LOOSENESS 
ENERGY VALUES (V^2) 





Frequency Good Damaged 
Shaft Spectrum 41.0 28.0 
Shaft Frequency 1.07 0.83 
tst Harmonic 2.31 1.22 
2nd Harmonic 1.23 2.40 
Upper Spectrum 31.4 26.31 


As can be seen in the table above, the energy content approximately doubled at the second 
shaft harmonic. This is a very encouraging result for the purposes of using the Pseudo 
Wigner-Ville Distribution energy as a fault detection tool. 
2. Gear Tooth Damage 

This particular fault was imposed on the gear model by shaving a gear tooth on 
the driver gear. The degree of damage was much easier to predict in this case and made this 
particular fault an easier one to control. The gear tooth that was damaged had approximately, 
1/8'^ of the leading edge removed. Shown in the table below are again the results of 
averaging two good runs and two post-damage runs and comparing the respective energy 


values at the different frequencies. For this case, just the frequencies associated with the gear 
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mesh are presented, since it is in these monitored frequencies that a fault involving the 


damage of a gear train would present itself. 


Table II. GEAR TOOTH DAMAGE 





GEAR TOOTH DAMAGE 
ENERGY VALUES (V^2) 


MESH FREQUENCY GOOD DAMAGED 


Fundamental 17.14 154.4 
1st Harmonic 0.53 3.44 
2nd Harmonic 0.16 0.48 


As shown in the above table, the energy content greatly increased at the fundamental gear 
mesh frequency. Additionaly, there was also increases at the harmonics and shown in Figures 
31 and 32 is a graphical presentation of the energy increase that occured due to the damaged 


tooth. 
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Figure 31. Gear Mesh Frequency (Pre-Damage) 
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GEAR TOOTH DAMAGE 
GEAR MESH INPUT ONE 
FILTERED 1250 - 1450 HZ 





Figure 32. Gear Mesh Frequency (Post-Damage) 
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This chapter has shown two examples of machinery damage and the resulting changes 
of the Pseudo Wigner-Ville Distribution energy. By applying the window and energy 
requirements discussed in Chapters III and IV, the Wigner- Ville energy has yielded a definite 
change to allow a means for machinery monitoring. The Wigner- Ville energy will make a 
stable input to a Neural Network and provide a quick efficient method by pre-processing the 
machinery signal with the Wigner- Ville Distribution as shown. This pre-processing will allow 
the number of inputs to the Neural Network to remain at a reasonable level and provide for 


a manageable, quicker hardware system. 


64 


VI. CONCLUSIONS 


The Pseudo Wigner-Ville Distribution has been shown to provide an excellent 
characterization of an input signal. Based upon the research conducted, it will also provide 
a means for machinery condition monitoring and diagnostics. The research has also shown 
that in order to use the Wigner-Ville Distribution in a machinery monitoring scheme, the 
following conditions or conclusions must be met and understood. 


e The smoothing window size applied must remain constant for purposes of comparing 
various signals. 


e The individual time and frequency lengths of the smoothing window may be adjusted 
to provide a better frequency resolution or constant amplitude. 


e The energy contained within the Pseudo Wigner-Ville Distribution of a signal is 
proportionaly time dependent. 


* The energy of the Distribution is proportional to the mean square value of the input 
signal and is not adversely affected by the cross-terms. 


e Lastly, the damaged gear model conditions showed a definite change in the energy of 
the Pseudo Wigner- Ville Distribution at the respective monitoring frequencies. 
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УП. RECOMMENDATIONS FOR FUTURE RESEARCH 


The Pseudo Wigner- Ville Distribution computer program that was established here at 
the Naval Postgraduate School is very stable. The pronounced peaks at the transition points 
are now definable and additionally, with the energy of the Pseudo Wigner- Ville Distribution 
now being a definable and understood quantity, the next additional research in this area may 
include: 


e Establish a baseline or data bank of "good" Wigner- Ville Distribution data for the 
present gear model. 


e Conduct additional runs of various levels of damage of gear model components for 
comparison with good data. 


e Establish a Torpedo Ejection Pump baseline of good data as with the gear model. 


e Continue with the research of linking the Pseudo Wigner- Ville Distribution as an input 
to a Neural Network Application. 
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APPENDIX COMPUTER PROGRAMS 
The following programs have all been developed here at the Naval Postgraduate School 
and used in the conduct of the research presented. The language used is FORTRAN 77 with 


the plotting programs using calls to CA-DISSPLA software. 


A. WVD3 COMPUTER PROGRAM 


PROGRAM WVD3 


* 
LE E EE EE EE 5555555 5 5 55 55 555 5 55 5 5 5 5 5 5 5 5 5 5 3 5 3 5 5 5 5 5 5 5 5 35 e fe e Жс fe ee d d dx 
* 


* THE FOLLOWING INDIVIDUALS HAVE BEEN INSTRUMENTAL IN THE 
* DEVELOPMENT OF THIS PROGRAM AND PREVIOUS VERSIONS: 


* 


Graham Rossano 
Jose G. Calahorrano 
Scott G. Spooner 


Prof. Y.S. Shin 


етте ез т ҰҚ + HH OF 


+ + њ»+ »+ + vg » * * 


De fe fe e e e e e De e de e e ee ee e De fee ee ee fe ee ce ee e ee ee ee fe Ve e Ve ede e fee de ec eh de d de x 


This program calculates the Wigner Distribution 
Function of a signal. It includes the option of applying a 
smoothing function and reducing the distribution to different 
sizes while writing the results to files for subsequent 
plotting. 


+ + ў ў д ях ж 


De fe e fe fe e e e fe e e К e e e ee efe ee e ee e e e De fe e e e ee ee ee ee ee ee e e ee e e ee de we ke we e ke e 


* * 
* VERSION 3 * 
* * 
LE E E E E 5 5 55 55 55 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 5 5 35 5 5 5 5 5 5 5 5 Z 3 5 5 5 5 3% 5 5 5 5 5 3 5 5 5 2 9 
* * 
* VARIABLE DECLARATION * 
* * 


fe c d e e fe e che e fe e ee e nhe efe ee ee fe e e e ee efe de ce e e ee e ee e Ve de fe ch c chc e e e eee e de e de ee de e 
* 
55553433343 ЖМ setting number of data points œ œ a Tr Ve e e e Yc We We de de de Œe œœ 
* 
* 
* 


integer dp,dp2,mvopt,redopt,dimt,dimf ,mn,nn,rf,rdp,rdpe, 

: nf,mt,nf2,mt?2 

PARAMETER (DP = 512) 

real pi,tin(dp),ain(dp), rawa(dp), rwdf(dp*2,dp) ,rswdf (dp*2,dp), 
rawt (dp) ,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb,val,sval,wdf(dp*2,dp), coef ,df,AREA 

complex s(dp*2),dum,c(dp*5) , dum3 , dum2 

character*25 inname 


common/in/ mvopt, redopt ,dimt,dimf,mm,nn,rf,rdp,rdp2, 
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: nf,mt,nf2,mt?2 
common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 
: rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
: const,t,sum,sumb,val,sval,wdf,coef,df,area 
common/comp/ s,dum,c,dum3,dum2 
арг - dp*? 
ж 


LE E 3 8 5 8 3 5 5 5 5 5 5 5 5 3 3 3 3 5 5 5 5 5 5 5 8 5 3 5 5 h 5 5 5 5 3 h h h 5 5 5 3 5 h 5 5 2 3 5 5 5 КЕ 9 


* INITIALIZING VARIABLES * 


fe e de de de e de e ede dede de Me dede ee eee ee ee ee ee e e ee ee ee ee e ee e e dede ee dee dede dede 


5 --- Print description of program --- 
print* 
print*,' Program to calculate the' 
print*,' Pseudo Wigner-Ville Distribution' 
print* 


S --- Set initial values --- 


print*,' Enter name of signal input file' 
read(5,904) inname 
print* 


print*,' Do you wish to remove the mean value?' 
print*,' Enter 1 for yes - O for по! 
read(5,910) mvopt 

print* 


print*,' Input the reduction size desired' 
print*,' input 1 for 64 by 32 - 
print*,' input 2 for 128 by 64 ! 
print*,' input 3 for 128 by 128 : 
print*,' input 4 for 256 by 128 | 
read(5,910) redopt 

print* 


print*,' Input the smoothing window size desired! 
print*,' input the size for the frequency axis' 
read(5,910) nf 

print*,' input the size for the time axis' 
read(5,910) mt 


* 
* 
* 
t e fe e e ne ne e e ne c e nee e ne eic ne e e ne de e e e ene e n nee ec nee e en e Òe Òe de de ne ee ce fe e de e de dede d x 


* * 
* calculation of some common used constants * 
* * 


LE E E E d E 5 8 53 3535454535338543 3 35 5 5 d 8 3 3 3 35 33 35 3 3 5 3 3 3 5 4 3 5 3 35 3 5 3 3 5 5 3 3 5 3 3 2 4 


pi = 4.0 * atan(1.0) 


ІЧЫ И ӨК Жї Ө e Me Ye Ye Ye Ae Me Ye de Me Ye de de e Me Me Ye Me КККК КК ККК КККК КК 


* THIS IS THE CALCULATION PART OF THE PROGRAM - 


c e e e e de e ede e ede de ce de e de e c e e f de e ne ee fe ede dede dede ede e ede de dede de dede dede de dede de fe dde dede dex 


open(4, file=inname,status='old') 
rewind 4 
* PRELIMINARY CALCULATIONS: 
CALL DATAIN3 
write(6,*) , ' done with datain3' 


CALL DTCALC3 
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write(6,*) , ' done with dtcalc3' 

CALL MEAN3 

write(6,*) , ' done with mean3! 

write(6,*) , ' mean value removed =', meanv 


* SIGNAL MOOIFICATIONS: 


* Window application (modified hamming window) 


CALL HAMMG3 
write(6,*) , ' done with hammg3'! 

x Conversion of real signal to an analytical one 
CALL ANAL3 


write(6,*) , ' done with апа(3' 


* CALCULATION OF THE WIGNER VILLE DISTRIBUTION: 
CALL WIGH3 


write(6,*) , ' done with wigh3' 


* REDUCTION & SMOOTHING OF THE WIGNER VILLE DISTRIBUTION: 
CALL REDUCE3 
write(6,*) , ' done with reduce3' 
CALL SMOOTH3 


write(6,*) , ' done with smooth3' 


їх й к ко ко ко о о ко о ік ік ік ік о de ко de ко ік ік de de de ко e ee e e ce ee te een ede e de de fe de ke de f k k k n e 


* * 
5 WRITING OF THE DIFFERENT ARRAYS GENERATED f 
x TO FILES FOR PLOTTING PURPOSES. * 
* * 


Tc Ïc he e he e e ee e nene eee eee een eee eee ee ee ke ke e edd d d d x x 

* 

* OPENING OF THE OUTPUT FILES 

* 
OPEN(UNIT=7, FILE="RAWTIME.OUT!, STATUS='NEW' ) 
OPEN(UNIT=8, FILE='"WNDWTIM.OUT', STATUS='NEW!) 
OPEN(UNIT=9, FILE='RWDF.OUT', STATUS='NEW' ) 
OPEN(UNIT=10, FILE='RSWDF.OUT' , STATUS='NEW') 
OPEN(UNIT=11, FILE="WOF.OUT', STATUS='NEW') 


* WRITING OF RAW AND WINDOWED TIME SIGNAL TO FILES 
* 


DO 300 І = 1,0Р 
WRITE(7,930) ТІМ(І),ВАМА(1) 
WRITE(8,930) ТІМ(І),АІМ(І) 
300 CONTINUE 
* 


* WRITING OF REDUCED/UNSMOOTHED WVD TO FILE 
* 


DO 400 1 = 1,DP/mm 


DO 400 J = 1,DP2/nn 
WRITE(9,906) RWOF(J,I) 
400 CONTINUE 
* 


* WRITING OF REDUCED & SMOOTHED WVD TO FILE 
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DO 500 I s 1,DP/mm 
DO 500 J = 1,DP2/nn 
WRITE( 10,906) RSWOF(J,1) 
500 CONTINUE 
* 


* WRITING OF ENTIRE WOF TO A FILE 
* 
00 600 1 = 1,DP 
DO 600 J = 1,0Р2 
WRITE(11,906) WOF(J,1) 
600 CONT I NUE 


* Format statements 


904 format(a25) 

906 format(2X,916.8) 

910 Тогта (16) 

950 format(2x,916.8,5x,916.8) 


END 
* SUBROUTINES 


include 'ANAL3.INC' 
include 'CORRS.INC' 
include 'DATAINS.INC' 
include 'DTCALC3.INC' 
include 'FFT3.INC' 
include 'HAMMG3.INC' 
include 'MEANS.INC' 
include 'SMOOTH3.INC' 
include 'REDUCE3.INC' 
include 'WIGH3.INC' 
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B. SUBROUTINES OF WVD3 


Listed below in alphabetic order are the subroutines used in the main program WVD3. 


SUBROUTINE ANAL3 


* 
їх о ко о ко ко о ік k о о ко о о о о ко о о о о ко о о по ко ко то кпк о то то о dr r dr о зо бо r о r k k e k k ок ко ко 


* * 
* This subroutine converts a real signal to an * 
* analytic one * 
* * 


LÈ E E 2 8 2 R 8 B B E о о о Мо о Мо о E X E 8 E 5 5 5 5 5 2 à à 8 8 8 E 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 à à ë d 2 r 
* 
* 
* 
integer dp,dp2,mvopt, redopt,dimt,dimf,mm,nn,rf,rdp,rdpe, 
: nf,mt,nf2,mt2 
PARAMETER (DP = 512) 
real pi,tin(dp),ain(dp),rawa(dp) , rwdf(dp*2,dp), rswdf (dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb, val, sval,wdf(dp*2,dp), coef ,df,area 
complex s(dp*2),dum,c(dp*5 ) , dum3 ,dum2 
character*25  inname 


common/in/ mvopt,redopt,dimt,dimf,mm,nn,rf ,rdp,rdp?2, 
: nf,mt,nf2,mt2 

common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 

: rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
: const, t,sum,sumb,val,sval,wdf,coef,df,area 
common/comp/ s,dum,c,dum3 ,dum2 

арг = ар*2 


до 100 1=1,ар 
sum-0.0 


do 200 j=1,dp 
sumb=0 .0 


if(i-j.eq.0) go to 200 
nzi-j 
val=pi*n/2. 
5уа| з5іп(ма| ) 
sumb-ain( j)*sval*sval/val 
200 sum-sum* sumb 
s(i)semplx(ain(i),sum) 


100 continue 


END 


7] 


SUBROUTINE DATAIN3 


T fede de ede hue ede re ede he ede e ehe eee dee eee ee e e cheer e he efe de de de dede eh d dœ 
* 


This subroutine assumes the input file is in the following * 
format: time amplitude (2x,e16.8,5x,e16.8) with an end of * 
file indicator of a last entry 9999,,9999, 5 


* 
"Üe fe he e de de ee ee e eee ede КККК КККК КККК ККК ККК КККК ККК ecce fe dede ded ЖКК 
* 
* 
* 


+ » * я я 


integer dp,dp2,mvopt ,redopt ,dimt,dimf,mm,nn,rf,rdp,rdp2, 

: nf,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi,tin(dp),ain(dp), rawa(dp), rwdf(dp*2,dp), rswdf(dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del 1,del2, 
const,t,sum, sumb, val , sval ,wdf(dp*2,dp),coef , df , AREA 

complex s(dp*2) ,dum,c(dp*5 ) , dum3 , dum2 

character*25 inname 


common/in/ mvopt ,redopt ,dimt ,dimf,mmn,nn,rf,rdp,rdp2, 
: nf,mt,nf2,mt2 

common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 

: rawt,dtsum,delt,dt,asum,meanv,mtime,del 1,del2, 
: const,t,sum,sumb, val, sval ,wdf,coef,df,area 
common/comp/ s,dum,c,dum3, dum2 

dpe = dp*2 


* 


********cimple Loop to read in time & amplitude*****s*sxswwwx 
do 100 j=1,dp 
read(4,902) rawt(j) , rawa(j) 
tin(j) = rawt(j) 
аіп()) = rawa(j) 
100 cont inue 


е FORMAT STATEMENT 
902 FORMAT (2X ,E16.8,5X,E16.8) 


200 END 
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SUBROUTINE DTCALC3 


* 
she e e he eee e ee e ee refe hee ree re ecce cc e rec c e e x 


* * 
* This subroutine calculates the delta t of the signal * 
* * 


а e e he he he ee eee ee c eere e el ce e А С Б о кс о e e x 

* 

* 

* 
integer dp,dp2,mvopt,redopt,dimt,dimf,mn,nn,rf,rdp,rdp2, 
: nf ,mt,nf2,mt2 
PARAMETER (DP = 512) 

real pi,tin(dp),ain(dp),rawa(dp),rwdf(dp*2,dp) ,rswdf(dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb, val, sval ,wdf(dp*2,dp),coef , df , AREA 

complex s(dp*2),dum,c(dp*5) , dum3 , Фито 

character*25  inname 


common/in/ mvopt, redopt,dimt,dimf,mm,nn,rf,rdp,rdp2, 
: nf ,mt,nf2,mt2 
common/rl/ pi,tin,ain, rawa, rwdf ,rswdf, 
: rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del?, 
: const,t,sum,sumb,val,sval,wdf,coef,df,area 
common/comp/ s,dum,c,dum3,dum2 
dp2 = dp*2 

* 


* 
ic she he e hc eee e ee eere ee eene ко о о о с о ко e о ік оо оо ік e de о о о e x 


dtsum = 0.0 
до 100 1 = 1 , др-1 


delt = tin(i*1) - tin(i) 
dtsum = dtsum * delt 


100 continue 


dt = dtsum / (dp-1.) 


END 
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SUBROUTINE FFT3 


* 
hr she їйї n n e e n n йа ғай n n n d n n n x Жї 


* * 
* This subroutine calculates the Fast Fourier Transform * 
* * 


hr vir e e ie e e rr e nr n rr err dde e e de e rrr er d d e rr d e dn n dn d dn v 
* 


* 
* 
* 


Безе dp, dp2 ,mvopt, redopt, dimt, dimf,mm,nn,rf,rdp,rdp2, 
nf ,mt,nf2,mt2 
' PARAMETER (DP - 512) 
real pi,tin(dp),ain(dp),rawa(dp),rwdf(dp*2,dp),rswdf (dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del 1,del2, 
: const,t,sum,sumb, val ,sval , df (dp*2,dp) ,coef , df 
complex s(dp*2),dum,c(dp*5) ,dum3, dum2 
character*25  inname 


common/in/ mvopt,redopt,dimt,dimf ,mn,nn,rf ,rdp,rdp2, 

: nf ,mt,nf2,mt2 

common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 

Е rawt ,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb,val,sval ,wdf,coef,df,area 

common/comp/ s,dum,c,dum3,dum2 

dp2 = dp*2.0 


+ + 


const=dp2 
val=alog(const)/alog(2.)+.1 


j=1 
do 40 i=1,dp2-1 
if (i.ge.j) go to 10 
dum3=c ( і) 
c(j)=c(1) 
c ( ï )=dum3 
10 k=dp 
20 if (k.ge.j) go to 30 
j=j-k 
k=k/2 
go to 20 
30 j=j+k 
40 continue 
do 70 n=1,val 
coe f=2**n 
coef=coef/2 
dum2-cmplx(1.,0.) 
dum-cmpl x(cos (pi /coef ), -sin(pi/coef)) 


do 60 j=1,coef 
do 50 izj,dp2,2*coef 
1i=i+coef 
dum3=c (11 )*dum2 
c(ii)=c(i)-dum3 
c(i )=c(i)+dum3 
50 cont inue 
dum2=dum2* dum 
60 cont inue 
70 | continue 


END 
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* 


SUBROUTINE HAMMG3 


ir she cr he й о ок ко її її її її її її її їй її r r ir ar ar air kr air ir kr ar dr ar akr kr ke r k r r t t r rr 


* 


* 


* This subroutine applies a modified hamming window  * 


* 
* 


to the signal ain(t) * 


* 


ft fr he іо he ік he he о hr rr nre nr о ік о e de e ко e e ко e ік ХК 


* 
* 


100 


integer dp,dp2,mvopt, redopt, dimt,dimf ,mm,nn,rf,rdpo,rdp2, 


nf ,mt,nf2,mt2 


PARAMETER (DP = 512) 


real pi,tin(dp),ain(dp), rawa(dp) , rwdf (dp*2,dp), rswdf (dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb,val,sval ,wdf(dp*2,dp), coef, df,area 


complex s(dp*2),dum,c(dp*5) , dum3 , dum2 


character*25 inname 


common/in/ mvopt,redopt,dimt,dimf,mm,nn,rf,rdp,rdp2, 
nf,mt,nf2,mt2 


common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 


rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del?2, 
const,t,sum,sumb, val, sval ,wdf,coef,df, area 
common/comp/ s,dum,c,dum3,dum2 
dp2 = dp*2 


mt ime=dp*dt 
del 120. 1*mt ime 
del220.9*mt ime 
const=pi/del1 
do 100 j = 1, dp 
t = (j-1) * dt 
if (t.le.del 1) then 
ain(j) es ain(j) * (.5*(1.-cos(const*t))) 
elseif ((t.ge.del2).and.(t.le.mtime)) then 
ain(j)2ain(j)*(.5*(1.-cos(const*(mtime-t)))) 
else 
endif 
continue 


END 
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* 


SUBROUTINE MEAN3 


какаа НИКИ 


* 


* 


* This subroutine calculates and removes the mean value * 
* of the signal. x 


* 


* 


fe fe de de e ne e e le e e le e e e e ee e rh re e le le e le le e e ee e e e e e e e e de err d Me Me Me d x 


* 
* 


ь 


100 


200 


integer dp,dp2,mvopt, redopt, dimt, dimf ,mm,nn, rf, rdp, rdpe, 


nf ,mt,nf2,mt2 


“PARAMETER (DP = 512) 


real pi,tin(dp),ain(dp), rawa(dp), rwdf (do*2, dp), rswdf (dp*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb, val, sval, wif (dp*2,dp), coef ,df,AREA 


complex s(dop*2),dum, c(dp*5) , dum5 , dum2 


character*25  inname 


common/in/ mvopt,redopt,dimt,dimf,mm,nn,rf,rdp,rdpe, 


nf,mt,nf2,mt2 


common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 


rawt,dtsum,delt,dt,asum,meanv,mtime, del 1, del2, 
const,t,sum,sumb, val ,sval,wdf,coef,df,area 
common/comp/ s,dum,c,dum3,dum2 
dpe = dp*? 


asum = 0.0 

do 100 i = 1, dp 
asum = asum + ain(i) 

continue 

meanv = asum / dp 

if (mvopt.eq.1) then 
do 200 1 = 1, dp 
ain(1) = ain(1) - теапу 
continue 

endif 


END 
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SUBROUTINE REDUCE3 


* 
he e hr he he she Ааа e e e n e he e e he e а а ee he he Ye dèr dr dir hee e e e e e e e e v Я d n 


* * 
* This subroutine reduces the wigner ville 5 
* to a workable size. * 
* * 


sir sr ік ік e e e ко ко ко о n nhe о he e e he іо e e rr e e e ehe ко ко he he he e e ко ко e ко he he e ХИ ЙЛЖЖЖЛЭИЖ 


* 
* 


possen do, dp2,mvopt ,redopt ,dimt,dimf ,mm,nn,rf,rdp,rdp2, 
nf ,mt ,nf2,mt2 
PARAMETER (ОР = 512) 
real pi,tin(dp),ain(dp),rawa(dp),rwdf(dp*2,dp),rswdf(dp*2,dp), 
: rawt(do),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
: const, t,sum,sumb, val, sval, wdf(do*2,do), coef ,df,area 
complex s(dp*2),dum,c(dp*5) , dum3, dum2 
character*25 inname 


common/in/ mvopt,redopt,dimt,dimf,mm,nn,rf ,rdp,rdp2, 
: nf,mt,nf2,mt2 
common/rl/ pi,tin,ain,rawa,rwdf,rswdf, 


: rawt,dtsum, delt, dt ,asum,meanv ,mtime,del1,del2, 
5 const, t, sum, sum, val ,sval ,wdf, coef,df, area 
common/comp/ s,dum,c,dum3, dum2 

dp2 = dp*2 


P d 


IF (REDOPT.EQ.1) THEN 


2 = 64 
Т = 32 
dk (REDOPT.EQ.2) THEN 
Е = 128 
E = 64 
ELSEIF (REDOPT.EQ.3) THEN 
КЕ = 128 
RT = 128 
ELSE 
RF = 256 
RT = 128 
ENDIF 
mn = дрг / КЕ 
mm = dp / RT 
[20 
do 100 j = 1, др, т 
l=l+1 
к= 0 
до 100 1 = 1, dp , mm 
k=k +1 


rwdf(l,k) = wdf(j,1) 
100 continue 


END 


pu 


SUBROUTINE SMOOTH3 


* 
e e fe e e e ко оо о оо e ко e he e e de dee eee e de ee e ee ede ee dee e e e ко e e e e eee eee eed ко ко 


* 


This subroutine reduces and smooths the WDF along both 
the time and frequency axes for purposes of plotting 
clarity. 


+ + + + 
+ + + + + 


ІІ 35333553 32233332 3 5 3 3 3 5 h 5 5 3 5 3 3 5 h 3 D 3 5 5 Z 2 5 3 5 5 2 8 5 5 5 5 2 2 2 
* 


real hg(-100:100, -100:100) 


integer dp,dp2,mvopt,redopt ,dimt,dimf ,mn,nn,rf,rdp,rdp2, 

3 nf ,mt,nf2,mt2 

PARAMETER (DP = 512) 

real pi,tin(dp),ain(dp), rawa(dp), rwdf (dp*2, dp), rswdf (dp*2,dp), 
rawt (dp) ,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum, sumb, val, sval,wdf(dp*2,dp) ,coef ,df, area 

complex s(dp*2),dum,c(dp*5) , dum3, dum2 

character*25  inname 


common/in/ mvopt,redopt,dimt,dimf,mn,nn,rf,rdp,rdp2, 

: nf ,mt,nf2,mt2 

common/ri/ pi,tin,ain, rawa,rwdf ,rswdf, 
rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del?2, 
const,t,sum,sumb, val ,sval,wdf,coef,df ,area 

common/comp/ s,dum,c,dum3 , dum2 

dpe = dp*2 


* 


dfz1./(4.*dp*dt) 


rdp-dp/mmn 
rdp2=dp2/nn 


nf2=nf*2 
mt2=mt*2 


valz1./((2.*pi)**2*nf*mt) 


do 20 j=-mt2,mt2 
do 10 i=-nf2,nf2 
coefz-((j*j)/(2.*mt*mt)*(i*1)/(2.*nf*nf)) 
hg(i, j )=val*exp(coef ) 
10 continue 
20 continue 
do 100 j=1,rdp 
do 100 iz1,rdp?2 
100 rswdf(i, j)=0.0 


* 
* 
ії з 0 
ро 500 N=1,DP2,nn 
jj = 0 
іі з 111 
DO 450 M = 1,DP,mm 
jj = jj*1 


DO 400 Kzm-MT,m*MT 
DO 350 L = n-nf, mnf 


IF (L.LT.1) THEN 


RSWDF(11,Jj) = RSWDF(11,Jjj) + O 
ELSEIF (K.LT.1) THEN 
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350 
400 


450 
200 


Ё5МОР(11,))) = 85м0Е(11,1))) 90 
ELSEIF (L.GT.DP2) THEN 
RSWOF(ii,jj) = RSWOFCii,jj) + O 
ELSEIF (K.GT.DP) THEN 
RSWOF(Ii,jj) = RSWOFCIii,jj) + 0 
ELSE 
RSWOF(ii,jj) = RSWOFCii, jj)+WOF(L,K)*HG(L-N,K-M) 
ENDIF 


CONT INUE 
CONT INUE 
К5МОҒ(11,)))-Е5МОҒ(11,))) 
CONTINUE 
CONTINUE 


END 
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SUBROUTINE WIGH3 


* 
dc e d dc de d de de e ee De e e ee e de e de e ccc e e e e e ec e e e e ncc e e e ee e d e e e e de e Wœ 


* * 
x This subroutine calculates the WDF of the signal % 
* * 


Tc dc c c de de de c e e de e e e de e e e de e e x Ve Ye Ve de dx Ye De de Dc c Ve Ve e e ee e nce e e e e de e e e e e e e de x 
* 
* 
* 
integer dp,dp2,mvopt,redopt,dimt,dimf ,mm,nn,rf,rdp,rdp2, 
: nf ,mt,nf2,mt2 
PARAMETER (DP = 512) 
real pi,tin(dp), ain(dp), rawa(dp), rwdf (dp*2,dp), rswdf (dpo*2,dp), 
rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
: const,t,sum,sumb, val ,sval , df (dp*2,dp),coef , df , area 
complex s(dp*2), dum, c(dp*5 ), dum3 , dum2 
character*25 inname 


common/in/ mvopt, redopt,dimt,dimf,mm,nn,rf,rdp,rdp2d, 
: nf ,mt,nf2,mt2 
common/rl/ pi, tin, ain, rawa, rwdf,rswdf, 
rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
const,t,sum,sumb, val ,sval ,wdf,coef,df,area 
common/comp/ s,dum,c,dum3,dum2 
dp2 = dp*2 
до 100 j = 1, dp 
яааааал авто correlation of the signal corr3 
coef = 2.0 * dt 
do 300 1 5 1 , dp*1 
if(j.ge.1) then 
dum = s(j-i+1) 
else 
dum = cmplx(0.,0.) 
endi f 
c(i) = coef * (s(j+i-1)*conjg(dum)) 
if(i.ne.j.and.i.ne.dp+1) then 
с(дрг-1+2) = сопіа(с(і)) 
endi f 
300 continue 
call fft3 
do 200 i=1,dp2 
wdf(i,j)=real(c(i)) 
200 continue 
100 continue 


END 
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C. PLOTTING PROGRAMS 
Listed below are the two computer programs that were used to produce the Figures in 
this work. 


1. 2-D Plotting Routine 


PROGRAM SIGNALPLT 


* 
fe e e e e e e e 8 5 8 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 3 5 5 5 5 5 5 5 5 4 5 5 5 5 5 5 5 5 5 АХ 5 4 5 5 5 5 d К 
* 


5 THIS PROGRAM USES THE GRAPHIC PACKAGE DISSPLA TO 
Я PLOT THE A 2-D FUNCTION. 


* 
RARARAKKKCKKKKKKKKKKKKKKKAKAKAKAKKAKKKKAKKKKKKKKKKKKKKEKEKEKEK 
* 
* — ***---DECLARING VARIABLES-- -*** 
* 

REAL XARAY(1000) , YARAY(1000) 

REAL X,Y,F,XSTART, XEND, YSTART, YEND, XSTEP 

INTEGER IPROPT,NPNTS, IX, 1Y 


***---[NITIALIZING VARIABLES AND SETTING GRAPH OPTIONS---*** 


» х к 


---set option for x starting & ending value--- 
WRITE(*,*),'DECIDE THE STARTING AND ENDING VALUES OF X' 
WRITE(*,*),'INPUT THE STARTING VALUE OF X' 
READ(*,*), XSTART 
WRITE(*,*),'INPUT THE ENDING VALUE OF X' 
READ(*,*), XEND 
2 ---set desired x step size--- 
WRITE(*,*), ‘DECIDE THE X STEP SIZE' 
WRITE(*,*),'INPUT THE X STEP SIZE! 
READ(*,*), XSTEP 

* ---set option for y starting & ending value--- 
WRITE(*,*),'DECIDE THE STARTING AND ENDING VALUES OF Y' 
WRITE(*,*),'INPUT THE STARTING VALUE OF Y' 
READ(*,*), YSTART 
WRITEC*,*), INPUT THE ENDING VALUE OF Y' 
READ(*,*), YEND 

E ---set option for number of points to plot--- 
WRITE(*,*),'DECIDE THE NUMBER OF POINTS TO PLOT’ 
WRITEC(*,*),' INPUT THE NUMBER OF POINTS! 
READ(*,*), NPNTS 

5 ---set option for line style (IMARK)--- 
WRITE(*,*),'DECIDE THE IMARK STYLE DESIRED' 
WRITE(*,*),* INPUT THE IMARK STYLE' 
READ(*,*), IMARK 

5 ---set option for grid style--- 

WRITE(*,*), ‘DECIDE THE GRID STYLE DESIRED' 

WRITE(*,*),'INPUT THE X GRID STYLE! 

READ(*,*), IX 

WRITE(C*,*5),'INPUT THE Y GRID STYLE' 

READ(*,*), IY 


* ---set option for output device--- 
WRITE(*,*),'YOU WILL NOW DECIDE WHERE TO SEND THE OUTPUT! 
WRITE(*,*), ENTER O FOR OUTPUT TO THE SCREEN! 
WRITE(*,*),'ENTER 1 FOR OUTPUT TO THE LASER' 
READ(*,*), IPROPT 
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CALL PDEV('LNO3', IEER) 
IF (ІРКОРТ.Е0.0) THEN 
CALL PGPX 
ELSE 
CALL LNO3I 
ENDIF 


* * * * 


**FUNCTION EXECUTION TO DETERMINE PLOTTING POINTS** 


OPEN (20,FILEZ'SINA00256A2.DAT' , STATUS 'OLD'! 
CALL F(XARAY , YARAY , NPNTS , XSTART , XEND ) 


+ ж 


* — ***---LEVEL ONE COMMANDS- - -*** 
CALL HWROT( 'COMIC') 
CALL PAGE(11.0,8.5) 
CALL NOBRDR 
CALL AREA2D(9.0,6.0) 


» 


***---LEVEL TWO COMMANDS- - -*** 


CALL SWISSB 

CALL HEADIN('400 HZ SINE МАМЕ5',100,1.5,2) 

CALL HEADIN('SIGNAL LENGTH = 0.128 SEC / AMP = 2.0$',100,1.5,2) 
CALL SWISSM 

CALL XNAME('TIME (sec)$!,100) 

CALL YNAME('AMPLITUDES' , 100) 

CALL GRAF(XSTART,XSTEP, XEND, YSTART, 'SCALE', YEND) 


» 


***---LEVEL THREE COMMANDS---*** 


CALL CURVE(XARAY, YARAY,NPNTS, IMARK) 
CALL FRAME 
CALL GRID(IX,IY) 


* 


***---CLOSING COMMANDS- - -*** 
CALL ENDPL(Q) 


CALL DONEPL 
END 


---FUNCTION SUBROUTINE TO DETERMINE PLOTTING POINTS--- 


* * * * 


SUBROUTINE F(XARAY,YARAY,NPNTS, XSTART , XEND) 


REAL XARAY (NPNTS) , YARAY (NPNTS) , XSTART,, XEND, STEP, XTEMP 
INTEGER NPNTS 
* ---READ IN THE XARAY AND YARAY VALUES--- 
DO 100 I = 1,NPNTS 
READ(20,*) ХАВАҮ(1), ҮАВАҮ(1) 
100 CONTINUE 
END 
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2. 3-D Plotting Routine 


PROGRAM WOFPLTVER3 


* 
fe e fe e he he de he ee eee ede e e he ecce e e dc e rece e ee ee eec e dede e de de d d n 
* 


5 THIS PROGRAM USES THE GRAPHIC PACKAGE DISSPLA TO 
е PLOT A 3-D FUNCTION. 


* 
*fe de de e e e eee e ee ee eee ee dee ee ede e ede dede e ede ede de dee dde de 3 8 6 n n n 8 
* 
* — ***---DECLARING VARIABLES-- -*** 
* 
REAL RWDF(32768) 
INTEGER IPROPT,IXPTS,IYPTS 
* 
* — ***---[NITIALIZING VARIABLES AND SETTING GRAPH OPTIONS---*** 
OPEN(15, FILE='"RSWOF.OUT',STATUS='OLD') 


DO 100 I = 1,32768 


READ(15,*) RWDF(I) 
100 CONTINUE 
* 


call pdev('ln03', ieer) 
**FUNCTION EXECUTION TO DETERMINE PLOT** 


E**---LEVELSONE: COMMANDS- --*** 


+ + + ++ 


CALL PAGE(8.5,11.0) 
CALL AREA2D(8.0,7.0) 


+ + 


***---LEVEL TWO COMMANDS---*** 


CALL VOLM3D(6.,6.,4.) 
CALL HEADIN('200 & 600 Hz SINE WAVES', 100, -1.5,2) 
CALL HEADIN(' FREQUENCY RESOLUTION ; 5 X 35$',100,-1.5,2) 
CALL X3NAME ( ' FREQUENCYS' , 100) 
CALL Y3NAME(' ',1) 
CALL 23NAME( 'AMPLITUDES' , 100) 
CALL VUANGL(-90.,0.,30.) 
CALL GRAF3D(0.,1250./5.,1250.0, 
: 0.,.2048/2.,0.2048, 
0., 'SCALE',.01) 


ха ***---L EVEL THREE COMMANDS---*** 


CALL SURMAT(RWOF, 1,256, 1,128,0) 
CALL FRAME 


*  #**---CLOSING COMMANDS---*** 


CALL ENDPL(O) 
CALL DONEPL 


END 
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